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An elementary statistical observation identifies generalizations of the Fuoss approximation for the 
probability distribution function that describes ion clustering in electrolyte solutions. The simplest 
generalization, equivalent to a Poisson distribution model for inner-shell occupancy, exploits measur- 
able inter-ionic correlation functions, and is correct at the closest pair distances whether primitive 
electrolyte solutions models or molecularly detailed models are considered, and for low electrolyte 
concentrations in all cases. With detailed models these generalizations includes non-ionic inter- 
actions and solvation effects. These generalizations are relevant for computational analysis of bi- 
molecular reactive processes in solution. Comparisons with direct numerical simulation results show 
that the simplest generalization is accurate for a slightly supersaturated solution of tetraethylam- 
monium tetrafluoroborate in propylene carbonate ([tea] [BF4]/PC), and also for a primitive model 
associated with the [tea] [BF4]/PC results. For [tea] [BF4]/PC, the atomically detailed results iden- 
tify solvent-separated nearest-neighbor ion-pairs. This generalization is examined also for the ionic 
liquid l-butyl-3-methylimidazolium tetrafluoroborate ([bmim] [BF4]) where the simplest implemen- 
tation is less accurate. In this more challenging situation an augmented maximum entropy procedure 
is satisfactory, and explains the more varied near-neighbor distributions observed in that case. 



I. INTRODUCTION 



a-ion is empty of 7-ions can be obtained from 



Ion clustering has long been an essential ingredient 
of our physical understanding of electrolyte solutions 
at elevated concentrations?"^ To describe pairing of a 
counter-ion of type 7 with an ion of type a, we focus on 
the radial distribution of the closest 7-ion to a distin- 
guished a-ion. We denote that normalized radial distri- 
bution by g^)[r). A famous discussion of Fuoss^ arrived 
at the approximation 
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for a primitive model of 1-1 electrolyte as in FIG. 1. Here 
q is the magnitude of the formal ionic charges, cL Q is the 
distance of closest approach, e is the solution dielectric 
constant, 2p is number density of ions, and (fc/3)" 1 = T 
is the temperature. We propose and test generalizations 
of Eq. ([lj in the following. 

Several complications of the distributions of near- 
neighbor ion-pairs motivate the generalizations that we 
develop. Firstly, ion-clustering can be particularly sen- 
sitive to non-ionic interactions. Com pariso n (FIG. 1) of 
atomically-detailed simulation resulta^ES with those of 
a corresponding primitive modeP^ straightforwardly ex- 
emplifies that point. Eq. ([I]) only treats classic ionic in- 
teractions. Secondly, even for primitive models the Fuoss 
approximation can be unsatisfactory (FIG. 2). Thirdly, 
nearest-neighbor distributions generally depend on which 
ion of an ion-pair is regarded as the central ion (FIG. 3). 
The radial distribution of the anion nearest to a cation is 
different from the radial distribution of the cation near- 
est to an anion, g^) (r) ^ g^ a (r). The approximation 

Eq. |T| is symmetric fl£^(r) = g^Jr) 

We are lead then to generalizations by recalling that 
the probability that a ball of radius r centered on an 
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P*r\ a (n = 0) = 47rp 7 J g { ^ a {x)x 2 dx 



(2) 



the assessment of the probability that the nearest 7-ion 
is further away than r. The simple estimate 



p 7 | Q (n = 0) « exp [- (n 7 | Q (r))) 



(3) 



with (rt 7 | Q (r)) = 47rp 7 fl g~ fa {x)x 2 dx, p 1 the density of 
7 ions, and g ia (x) the conventional radial distribution 
function, follows from the assumption of the Poisson dis- 
tribution for that probability. Evaluating the derivative 
of Eq. Q using Eq. ^ gives 



9y\l(r) ~ g ai {r) exp 



-47rp 7 / g ja (x)x dx 
'0 



(4) 



For g ai (r) = 1 (no correlations), this is th e Hertz dis- 
tribution that is correct for that caseP^^ We recover 
the Fuoss approximation with In g a ~ ( (r) ss — /3<7 a q 7 /er = 
/3q 2 /er for r > d Q7 , and zero (0) otherwise. This deriva- 
tion of the Fuoss approximation Eq. seems not to 
have been given before. Nevertheless, the suggested ap- 
proximation Eq. Q is a standard idea in the context of 
scaled-particle theories of the hard-sphere fluidP^ 

As discussed below, the Poisson result Eq. ^ follows 
from a maximum entropy development when the infor- 
mation supplied is the expected occupancy of the inner- 
shell EHED That information is sufficient if the occupancy 
n(r) is always low, i.e., rarely larger than one. Thus, 
in contrast to the Fuoss approximation, Eq. Q is cor- 
rect for small r because the expected coordination num- 
ber tends to zero then. For the same reason, the Pois- 
son approximation Eq. ([3| is correct at low electrolyte 
concentration, and even when the solvent is treated at 
atomic resolution. Furthermore, it is natural to guess 
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FIG. 1: Upper-panel: Ion-ion radial distribution functions for 
atomically-detailed simulatio n 112 * 13 ^ of tetraethylammonium 
tetrafiuoroborate in propylene carbonate ([tea][BF4]/PC) at 
T — 300 K, p — 1 atm, and the slightly supersaturated concen- 
tration of 1 mol/dm 3 . g^ N (r) is the radial distribution of the 
nearest B-neighbor of an N-atom. Lower panel: Results for 
a corresponding primitive model with dielectric constant and 
with ion charges and sizes matched to the [tea] [BF4]/PC case 
above. Specifically the model dielectric constant is e = 60, 

and d++ = 0.6668 nm, d = 0.6543 nm, d |_ = 0.45 nm 

are distances of closest approach for the hard spherical ions. 
The lower panel was produced by Monte Carlo simulation of 
a neutral system of 2 x 500 hard spherical ions in conventional 
cubical periodic boundary conditions at the same temperature 
and concentration as the results above, utilizing the Towhee^ 
package adapted to the present system. 



cation-anion chain or ring structures when ionic inter- 
actions drive well developed clustering. FIG. 1 shows a 
mean coordination number of less than two for counter- 
ion neighbors closer than about 0.5 nm, and supports the 
chain/ring picture of ion clusters formed. It is plausible 
therefore that a choice of inner-shell radii leading to small 
coordination numbers should validly describe important 
features of well-developed ion-clustering. 

For computational analysis of reactive bi-molecular en- 
counters in solution, identi ficati on of geometries of closest 
molecular pairs is critical p212U Because it is correct for 
low concentration and for small r in any case, Eq. Q 
should be regarded as the general resolution of those 
questions. 

When coordination numbers exceed one with reason- 
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FIG. 2: Comparison of the Fuoss approximation for g^J_(r) 
to Monte Carlo results of FIG. 1. 




FIG. 3: An example showing that the distribution of anions 
nearest to a cation is generally different from the distribution 
of the cations nearest to an anion. The dashed lines indicate 
the nearest distances of a cation to each of the three anions. 
The solid lines mark the nearest distances of a anion to each 
of the three cations. 

able probability, information on the expected number of 
pairs of counter-ions in the inner-shell should improve 
a maximum entropy model of these probabilities .UZHEl 
A maximum entropy model involving pair information 
would predict the g^L(t) 7^ 9^ a (f) asymmetry. For a 
1-1 electrolyte, the generalization Eq. Q is symmetrical 
in accord with the Fuoss approximation. The extent to 
which the observed asymmetry is significant gives an in- 
dication whether the Poisson approximation is adequate. 

In this work, the Poisson approximation (Eq. Q) is 
tested using three distinct simulation data sets. Two 
of these data sets have been noted already in consid- 
ering FIG. 1. Those calculations treated solutions of 
tetraethylammonium tetrafiuoroborate in propylene car- 
bonate, one at atomic resolution ([tea] [BF^/PC) and 
the other on the basis of a primitive electrolyte solution 
model over a range of concentrations. The third data 
set treated the ionic liquid l-butyl-3-methylimidazolium 
tetrafiuoroborate Qbmim] [BF 4 ]). To ensure the correct 
correspondence of the necessary simulation details with 
the results as they are discussed, those details are pro- 
vided in the captions of the figures providing the Simula- 
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tion results. 



II. RESULTS AND DISCUSSION 

For [tea] [BF^/PC, comparison (FIG. 4) of the numer- 
ical data with the approximation Eq. Q shows agree- 
ment over a distance range wider than the sizes of the 
molecules as judged by the radial distributions (FIG. 1). 
These near-neighbor distributions show bi-modal proba- 
bility densities with maxima at r rs 0.5 nm and 0.9 nm. 
These correspond, respectively, to a contact ion pair and 
to a solvent-separated near-neighbor ion-pair. Thus the 
Poisson approximation Eq. Q in this case includes solva- 
tion structure in characterizing inter-ionic neighborship. 
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FIG. 4: For the [tea] [BF 4 ]/PC case of FIG. 1, comparison 
of the numerical data with the approximation Eq. (|4j. The 
local maximum at r ~ 0.9 nm identifies solvent-separated 
nearest-neighbor ion-pairs. In this case, asymmetry of the 
two observed near-neighbor distributions is slight. That sug- 
gests that more elaborate maximum entropy models are un- 
necessary, and indeed the Poisson approximation is accurate. 
The embedded molecular graphic shows one of the solvent- 
separated nearest-neighbor BF4~...PC ...tea + structures 
observed. 

A plateau between r ss 0.5 nm and 0.9 nm in occupancy 
probabilities (FIG. 5) indicates saturation of counter-ion 
probability, and marks the inter-shell region. At the dis- 
tance r indicated by the vertical line, the coordination 
numbers n = 1, 2 predominate, supporting the idea of 
the formation of cation-anion chain and ring structures. 
The two sets of probabilities (FIG. 5) are qualitatively 
similar, reinforcing the symmetry of FIG. 4. 

Results (FIG. 6) for the primitive model of FIG. 1 ex- 
amine the sufficiency of the Poisson approximation over 
a broader concentration range for such models. The 
nearest-neighbor distributions are unimodal in this case. 
Correct at small r where the probability densities are 
highest and properly normalized, the Poisson approxi- 
mation Eq. Q is accurate over the whole range shown. 

Another example is the ionic liquid [bmim] [BF4] , with 
molecular structure shown in FIG. 7 and radial distri- 
bution functions in FIG. 8. The Poisson approxima- 
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FIG. 5: Occupancy probabilities as functions of the observa- 
tion sphere radii r for the [tea] [BF4]/PC case of FIG. 1. Upper 
panel: probabilities for occupancy by the B-atom of the BF4 - 
anion of the inner-sphere of the N-atom of tetraethylammo- 
nium cation (tea + ). Lower panel: probabilities for occupancy 
by the N-atom of the cation of the inner-sphere of the B-atom 
of the anion. The curves lowest in each panel are similar, 
showing symmetry displayed also in FIG. 4. 



tion (Eq. Q) agrees with the observed 3c2| B ( r ) at short 
range and displays a second maximum characterizing 
non-contact nearest neighbors, though in this case there 
is no additional solvent. The near- neighbor B|C2 dis- 
tribution (FIG. 9), on the other hand, lacks a second 

maximum. Thus <?c2| B ( r ) and 5B|C2( r ) ^ or i° n i c liquid 
[bmim][BF4] display the generally expected asymmetry. 
This asymmetry is also reflected in occupancy probability 
profiles (FIG. 10). More general theoretical models are 
required for such cases, and we return to that theoretical 
discussion now. 



III. MAXIMUM ENTROPY MODELING 



The Poisson distribution p (n) — 



-(n) 



jn\ de- 



scribes random occupancy consistent with the informa- 
tion (n) — (n(r)) . Considering the relative entropy, 



V({p(n)})= -^2p{n)\n 
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FIG. 6: Dashed curves are Monte Carlo results for the prim- 
itive model of FIG. 1 and the solid lines are the Poisson ap- 
proximation, Eq. Q. From top to bottom, the distinct cases 
correspond to concentrations 0.01, 0.05, 0.10, 0.20, 0.40, and 
0.80 mol/dm 3 , T = 300K for each case. For the highest con- 
centration, the system size is 2x400 ions. For all other cases, 
the system size is 2x200 ions. At the lowest concentration 
here the distribution of the nearest neighbor g_h(r) is close 
to the full radial distribution function (r). 
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FIG. 7: Structures and atom labeling of cation l-butyl-3- 
methylimidazolium (bmim + ) (top) and anion tetrafluorobo- 
rate (BF4 - ). On the cation molecule, white, turquoise and 
dark blue balls stand for hydrogen, carbon, and nitrogen 
atoms, respectively. On the anion, white and green balls stand 
for fluorine and boron atoms, respectively. 



FIG. 8: Radial distribution function of C2 atom on bmim + 
cation and boron on BF4 - anion obtained from MD simula- 
tion. Force field parameters and partial charge of the atoms 
were taken from Andrade et alW* The initial unit cell, with a 
dimension of 40 x 40 x 40 A is uniformly packed with 190 
ion pairs using PackmolP^ MD simulation was performed us- 
ing AMBER10 at constant pressure (1 bar) and temperature 
(298. 5K). The system was first minimized, followed by 0.2ns 
equilibrium at a time step of 0.2 fs, then 1.3 ns equilibrium 
at a time step of 2 fs. Radial distribution functions were 
extracted from a production run of 3.0 ns. 
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FIG. 9: Nearest-neighbor distributions for the ionic liquid 
case of FIG. 8. In this case with no solvent, the asymmetry 
of the two observed near-neighbor distributions is marked. 
An outer-sphere nearest-neighbor is exhibited in one case but 
not the other case where the Poisson approximation Eq. Q 
is inaccurate. 



the Poisson distribution is a maximum entropy distri- 
bution satisfying the specific expected occupancy. If we 
have more information, e.g., the binomial moments^H^D 



7 / / — V (n — iV.1l 

J 7 ' n>0 " ' 



(6) 



we can seek the distribution which maximizes rj({p( n )}) 
and satisfies the broader set of information. 

With the binomial moments (Eq. (6])), the Poisson dis- 
tribution is seen to be correct if realized values of n are 
rarely bigger than one (1). If n is never 2 or larger, bi- 



nomial moments j > 2 vanish. When j > 2 binomial 
moments are small, and that is consistent with Poisson 
prediction that they are zero. This underlies our ob- 
servation above the the Poisson model, p(0) ~ e~^ of 
Eq. pi), is correct for small A. 

Beyond the mean occupancy, the next level 
of information is the pair-correlation informa- 
tion (n(r) (n(r) — 1) /2), the expected number of 
pairs of counter-ions in the indicated inner-shell. 
Carrying-out the maximization for the case that 
pair information is available induces the model 
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FIG. 10: Occupancy probabilities as functions of the ob- 
servation sphere radii r for [bmim][BF4] (FIG. 8). Upper 
panel: probabilities for occupancy by the B-atom of the BF4 - 
anion of the inner-sphere of the mid-C-atom of l-butyl-3- 
methylimidazolium (bmim + ) cation (labeled as C2 in FIG. 7). 
Lower panel: probabilities for occupancy by the mid-C-atom 
of the cation of the inner-sphere of the B-atom of the anion. 
These results illustrate the possibility of a structural motif 
of ion clusters as chains and rings, i.e. at the r = 0.43 nm 
distance of the vertical dashed line probabilities of 1 and 2 
neighbors predominate. They also demonstrate asymmetry 
of the distributions of near-neighbor distances in their depen- 
dence on choice of the central ion, i.e. the distribution of 
nearest anions to a cation is different from the distribution 
of the nearest cations to an anion. Since p(3) is larger in the 
lower panel than in the upper panel, the BF4 - anion is more 
likely to be a three-way junction in this analysis than is the 
bmim + cation. 



p(n) oc exp [— Xin — \2?i(n — l)/2] /n\, where Ai, 
and A2 are Lagrange multipliers adjusted to reproduce 
the information (n) and (n(n— l)/2). Explicitly 
addressing the normalization of these probabilities leads 
to 



p(n) 



-XxTl— A 2 n(n— 1)/ 2 



p(0) involves only the denominator of Eq. ([7]), and 
can be considered a partition function sum over occu- 
pancy states with n-dependent interactions and interac- 
tion strengths adjusted to satisfy the available informa- 
tion. The information required (FIG. 11) for this aug- 
mented maximum-entropy model is only subtly different 
for the two cases. Nevertheless, the results (FIG. 12) 
agree nicely with the observed asymmetry. 




FIG. 11: Pair information required for the augmented max- 
imum entropy prediction of Eq. pj. Solid curve: the num- 
ber of C2 pairs occupying a sphere of radius r on a B atom. 
Dashed curve: the number of B pairs occupying a sphere of 
radius r on a C2 atom. 
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FIG. 12: Analysis of the asymmetric n = probabilities 
of FIG. 10. Solid curves: direct numerical simulation as 
in FIG. 10; Dashed curve: Poisson-based approximation, 
Eq. Q; Dotted curves: augmented maximum entropy model 
utilizing the first two binomial moments, Eq. j8j. The two- 
moment maximum entropy model is qualitatively reliable and 
therefore gives a satisfactory explanation of the observed 
asymmetry. 
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IV. CONCLUSION 



( 8 ) Results for both the [tea] [BF 4 ]/PC (FIG. 1) and the 

ionic liquid [bmim] [BF4] (FIG. 8) identify a natural 



() 



clustering radius where mean coordination numbers are 
near two. This suggests arrangements of the closest 
neighbors leading to a structural motif of cation-anion 
chains and rings. In contrast to the atomically detailed 
[tea] [BF4J/PC results, a corresponding primitive model 
(FIG. 1) does not display those clustering signatures 
(FIG. 6). A generalization (Eq. Q) of the Fuoss ion- 
pairing model was obtained by recognizing that the Pois- 
son distribution is correct when the mean coordination 
numbers are low. On the basis of measurable molecu- 
lar distribution functions, this generalization also estab- 
lishes the distribution of molecular nearest neighbors for 
computational analysis of bi-molecular reactive processes 
in solution. This Poisson-based model is accurate for 
the [tea] [BF^/PC results, both for the primitive model 



and the atomically detailed case. For [tea] [BF^/PC, the 
atomically detailed numerical results and the statistical 
model identify solvent-separated nearest-neighbor ion- 
pairs. Distributions of nearest-neighbor distances typi- 
cally depend on which ion of a pair is taken as the central 
ion, i.e., the distribution of anions nearest to a cation 
is different from the distribution of the cations nearest 
to an anion. The Poisson-based model is not asymmet- 
ric in that way. The numerical data for the ionic liquid 
[bmim] [BF4] prominently show the expected asymmetry. 
That asymmetry can be treated by a maximum entropy 
model based on the expected number of pairs of counter- 
ions occupying the inner-shell of the central ion, infor- 
mation extracted from the simulations. 
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